/*
 * 2-dimensional Ising model, with Markov chain in Monte Carlo method
 * 	H = - J \sum_{ij} s_i s_j - muB\sum_i s_i,
 * where J is the interaction strength between two atoms, negative/positive for ferromagnetic/anti-ferromagnetic material, k is the Boltzman constant, T is the temperature, we set J as unit.
 */

#include<iostream>
#include<stdlib.h>
#include<time.h> 
#include<fstream>
#include<cmath>
using namespace std;

void Ising2d(int length, int L, double kT, double muB, double & sz_ave, double & e_ave, double & e2_ave){

	long i,j,k, Sz, S_neighbour, count=0;
	double J=1, E, dE, sz_sum=0, e_sum=0, e2_sum=0;
	int *spin=new int [L*L];//spin of each atom, which can only be 1 or -1
	bool flag_flip;

	for(i=0;i<L*L;i++) //initially the lattice is magnetizied.
		spin[i]=1;

	srand((unsigned)time(0));//initialize rand() with a random seed, according to the real time when this code runs.

	Sz= L*L;//initial value of total spin
	E=-J*2*L*L - muB * Sz;//initial value of total energy of the system

	ofstream fp("IsingModelVisualization.gnu");
	fp<<"set term gif animate delay 10\n";
	fp<<"set output \"IsingModel.gif\"\n";
	fp<<"unset xtics\n";
	fp<<"unset ytics\n";
	fp<<"unset border\n";
	fp<<"set title \"kT = "<<kT<<"J, n=0, <s>=0, <E>=0\"\n";
	fp<<"set dummy i,j\n";
	fp<<"do for[i=0:"<<L-1<<"]{\n";
        fp<<"\t do for[j=0:"<<L-1<<"]{\n";
        fp<<"        set arrow from j+0.5,i+0.2 to j+0.5,i+0.8 size 3,0.8 filled lw 1 lc rgbcolor \"red\"\n";
	fp<<"\t }\n";
	fp<<"}\n";
	fp<<"plot \"rectangle.txt\" with line lc rgbcolor \"black\" notitle\n";

	for(j=0;j<length;j++){

		for(i=0;i<L*L;i++){
			S_neighbour = spin[ (i>(L-1)*L-1 ? i-(L-1)*L : i+L) ] + spin[ (i<L? i+(L-1)*L : i-L) ]
				+ spin[ (i%L ==0 ? i+L-1 : i-1) ] + spin[ ((i+1)%L==0 ? i-(L-1) : i+1) ];
			dE= 2 * J * S_neighbour * spin[i] + 2 * muB * spin[i];

			flag_flip=false;
			if(dE<0){
			       flag_flip=true;
			}
			else{
				if( (double)rand()/RAND_MAX < exp(-dE/kT) ){// rand()/RAND_MAX uniformly distributes in [0,1], so it has possibility P to be less than P
					flag_flip=true;
				}
			}
		
			if(flag_flip){
				E += dE;
				spin[i] *= -1;//flip the spin
				Sz += 2*spin[i];// total spin and total energy changes
			}
		}
		// after one round, collects data

		if(j>0.1*length){
			//cout<<"sz="<<(double)Sz/L/L<<"\t e="<<E/L/L<<endl;
			sz_sum += (double)Sz/(L*L);
			e_sum += E/(L*L);
			e2_sum += E*E/(L*L);

			count ++;
			if( count% (length/100)==0){
				fp<<"set title sprintf(\"kT="<<kT<<"J, n=%d, <s>=%f, <e> = %f\","<<count<<","<<sz_sum/count<<","<<e_sum/count<<")\n";
				for(i=0;i<L*L;i++){
					fp<<"unset arrow "<<(i+1)<<"\n";
					if(spin[i]>0)
						fp<<"set arrow from "<<(i%L)+0.5<<","<<i/L+0.2<<" to "<<(i%L)+0.5<<","<<i/L+0.8<<" size 3,0.8 filled lw 1 lc rgbcolor \"red\"\n";
					else
						fp<<"set arrow from "<<(i%L)+0.5<<","<<i/L+0.8<<" to "<<(i%L)+0.5<<","<<i/L+0.2<<" size 3,0.8 filled lw 1 lc rgbcolor \"blue\"\n";
				}
				fp<<"replot\n";
			}
		}
	}
	fp<<"set output\n";
	fp.close();
	sz_ave = sz_sum/count;
	e_ave = e_sum/count;
	e2_ave = e2_sum/count;

	delete [] spin;
	//double C=1/kT/kT*( E2_sample/num_sample - (E_sample*E_sample/num_sample/num_sample) );
}

int main(){

	int length = 1E4, L=16;
	double sz_ave, e_ave, e2_ave, muB=0;

	double kT; cout<<"kT="; cin>>kT;
	if(kT<=0){
		cout<<"kT must be greater than 0.\n";
		exit(1);
	}
	Ising2d(length, L, kT, muB, sz_ave, e_ave, e2_ave);

	return 0;
}
